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Abstract 

The set of 3D inviscid primitive equations for the atmosphere is dimensionally reduced by a 
Discontinuous Galerkin discretization in one horizontal direction. The resulting model is a 2D 
system of balance laws where with a source term depending on the layering procedure and the 
choice of coupling fluxes, which is established in terms of upwind considerations. The "2.5D" 
system is discretized via a WENO-TVD scheme based in a flux limiter centered approach. We 
study four tests cases related to atmospheric phenomena to analyze the physical validity of the 
model. 

Keywords: primitive equations, layering, discontinuous Galerkin, upwind flux, WENO-TVD 
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1. Introduction 

The so-called primitive equations used to model atmospheric and ocean dynamics are not 
well-posed for any reasonable set of boundary conditions. This important, and relatively old 
result, was obtained by Oliger and Sundstrom in d. The problem of formulating a well-posed 
set of "primitive" equations is of obvious importance both from the mathematical and numerical 
side, and has recently attracted substantial research activity. One example of such research was 
undertaken in O, where the authors have formulated a type of dimensionally-reduced equations, 
called 2.5D equations, where well-posedness was analyzed. However, this paper uses linearized 
primitive equations and the layering procedure was performed via a Continuous Galerkin ap- 
proach using orthogonal piecewise linear functions. This imposes some limitations regarding the 
linear character of the system and the number of elements to be considered in the expansion. 
What we want to communicate here is an alternative derivation of a 2.5D set of full nonlinear 
equations and its approximation by a high-order finite volume scheme. We consider the analysis 
of well-posedness for these equations to be a far too difficult task. Instead we perform numerical 
experiments with well-known test cases normally used to test dynamical cores of atmospheric 
models |[3l|4l. Good results from such experiments are a strong indication that our formulation 
of a 2.5D set of equations is well-posed. The dimensional reduction procedure that we propose 
is based on a Discontinuous Galerkin approach (DG) [Sj ; this class of schemes have proved to 
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be successful in the resolution 2D atmospheric models f3l. We use the same basic principles 
to divide the domain into 2D layers where the equations are locally approximated and coupled 
by suitable fluxes. This approach leads to a system of 2D balance laws, which is solved via a 
WENO-TVD scheme |6|, using a flux limiter centered approach |7 | for space discretization and 
a Runge-Kutta TVD scheme for time discretization. Since the primitive equations are not well- 
posed, one may wonder what approach is taken in the many atmospheric and oceanic models 
that exist. The answer is that a pragmatic, and hence not strictly rigorous, approach is taken. 
In practice this means that one imposes boundary conditions that are not boundary conditions 
in the ordinary sense. An example is the use of relaxation zones near the boundary, where the 
solution in the inner of the domain is relaxed towards a certain boundary value. While this is 
stable numerically (in most cases), it can hardly be said to be satisfactory in the physical sense. 
And not in the mathematical sense either. 

The remaining part of the paper is organized as follows: In section 2 we present the full 
3D governing equations in conservative form for a dry atmosphere, and its "semi-discretization" 
in a spatial dimension so the equations are formulated in 2.5D as explained above. The finite 
volume-based numerical scheme is presented in detail in section 3, and in section 4 we present a 
number of numerical experiments using both well-known tests for atmospheric model dynamical 
cores, as well as new tests. The results are discussed for each test, and we summarize the results 
along with directions for further work in section 5. 



2. A dimensionally-reduced system for atmospheric dynamics 

This section is concerned with the derivation of a dimensionally-reduced set of equations 
describing the atmosphere. We first present the original three-dimensional system of primitive 
equations which is then layered via a DG ansatz in order to obtain a 2.5D model. 

2.7. A three-dimensional system of balance laws 

Our starting point is the set of equations modelling the evolution in time of a dry atmosphere. 
We impose conservation of mass, momentum and energy, and consider eff'ects of rotation and 
gravity, together with neglecting friction; this leads to a set of 3D inviscid primitive equations 
cast in conservative form 



dtQ+ d,T + dyQ + d,^ = (1) 

where 
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The dynamics are expressed in terms of the density of the fluid p, the 3D velocity field (u, v, w), 
the pressure P and the potential temperature 0, which relates to the thermodynamic temperature 
T via 

The system is closed by the equation of state for an ideal gas, 

V = Coipey, Co = (5) 

^ 

The model parameters are: the Coriolis parameter / = 1 x 10""^ [^"^], the gravitational accel- 
eration g = 9.^l[ms~^], the atmospheric pressure at sea level = 10^ [Pa], the gas constant 
for dry air = 2Sl[JK~^kg~^], the specific heat of dry air at constant pressure and volume 
Cp = l004[JK~^kg~^], the specific heat of dry air at constant volume Cy = lll[JK~^kg~^] and 
its ratio y = Cp/Cy = 1.4. Moreover, by defining the Exner pressure 



the expression for the total energy of the system is given by, 

e = CyOn +]^(u^ + + w^) + gz. (7) 

The system ([T]) governs the time evolution of the conserved variables in Q by the interaction be- 
tween the physical fluxes T , Q and the source term S . A more detailed representation of the 
atmospheric dynamics could include an equation for conservation of moisture that will afl'ect the 
thermodynamics and hence the dynamics (expressed by the fluxes). Physical parameterizations 
in a weather model, like surface parameterizations and condensation processes, together with 
coordinate transformations, like the extensively used terrain-following coordinates will generate 
terms appearing in the source S . In the next subsection we present a dimensionally-reduced ver- 
sion of the original 3D set of equations where we characterize the behavior of the variables in one 
spatial dimension, in a way such that the action of its associated flux is transferred to the source 
term. Although this can be seen as a step in the numerical discretization of the original set of 3D 
equations, our study will be focused on the physical significance of the dimensionally-reduced 
(or 2.5D) set of equations as a system of balance laws on its own. 

2.2. Towards a 2.5D model 

In O, the authors obtained a 2.5D model from the linearized primitive equations via a Con- 
tinuous Galerkin expansion of the variables in the j-direction. It has to be emphasized that the 
derivation made considerable use of the linearity of the system, and that the basis functions were 
a set of orthogonal piecewise linear functions; removing any of those aspects leads to a model 
whose complexity is not clearly related those of the underlying physics. On the other hand, there 
is an obvious limitation related to the trade-off' between orthogonality and locality of the basis 
functions; the aforementioned article considered two elements with a basis of three orthogonal 
piecewise functions, which is just a slight modification of the classical hat functions. An incre- 
ment in the number of elements, together with the orthogonality condition, will produce a set 
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of basis function with increasing support, significantly diff'erent from the analogous set of hat 
functions. The approach that we present pays a cost by removing the continuity requirement on 
the solution (in the classical Galerkin framework), and switches to a DG formulation that allows, 
in a straightforward and tractable manner, to deal with nonlinearity and locality at the same time. 
We start by defining our set of equations over D = [0, U] x [0, U] x [0, L^] x [0, Tf]. Initial and 
boundary conditions will be specified later. We consider a uniform partition in the y-direction 
into A^3; elements, i.e. 

Y = [0, = U Yt, Yt = [j,--i/2, yi^i/2] yi^m - yi-m = Aj, V/. (8) 

Next, in every element we consider a local classical Galerkin ansatz by multiplying eq. ([T]) by a 
test function V and integrating with respect to y, leading to 



X 



V^ldiQ+ d,r + dyQ + d.'H - S}dy = 0. (9) 



-1/2 

Mapping every element into the canonical element (-1, 1) by means of 
and, as usual, performing integration by parts, leads to 

^ ^\T{d,Q+d,r + d,^ -s}d^ + v^0\\ = 9 d^. (1 1) 

We consider a local set of basis functions and expand Q restricted to F/: 

^, zj) = Yj mM)' (12) 

The choice of the set of basis functions is an open question and will depend on the application 
and the complexity of the computations. In this particular case, in order to be able to carry 
on in a tractable manner the calculations for a nonlinear model, it is important that we chose 
a set of orthogonal functions, and therefore selecting a set of orthogonal polynomials (as the 
Legendre polynomials for instance) will allow us to preserve a reasonable model but also to 
include high-order approximations. Throughout this article however, we will restrict ourselves 
to the most basic case in this framework, i.e., to consider Legendre polynomials up to degree 
0, which is nothing but to consider a piecewise constant approximation of the variables inside 
every element. In 1 8 |, the piecewise linear extension of this procedure was performed in a similar 
framework, and poses no additional difficulties except for the increase in the computational cost 
associated with the amount of local expansion coefficients. When we consider 2 elements in the 
j-direction, Yi = [0, 0.5 U] and Y2 = [0.5 U, U] and the first Legendre polynomial 0o = 1, the 
2.5D system of equations reads, 

dtQi^d,r(Qi)^d,9((Qi) = Ay{g(y = 0.5Ly)-g(y = 0)}^S(Qi), (13) 



dtQ2+d,r(Q2)+d,^(Q2) = Ay{g(y = Ly)-g(y = 0.5Ly)}+S(Q2y (14) 
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At this point we must close the system by giving computable expressions for the evaluation of 
Q at the boundaries and the intermediate value. We first address the problem of the intermediate 
value, and then we expand this procedure to the boundary fluxes. Following the same ideas of 
the finite volume framework, we seek to express the intermediate state as a function of adjacent 
cells, i.e. 

0{y = O.5Ly)^ G{QuQ2). (15) 

As it is extensively described in the literature, the two main choices for such procedure are 
centered (averaged) values, or states obtained via upwind considerations. The typical centered 
fluxes can be interpreted as particular averaging operators in space and time, which does not fit 
in a proper manner in our scheme, as time integration will be considered only for the resulting 
2.5D model; therefore we opt for upwind fluxes to determine the coupling at the interface. We 
recall that the flux Q in the Euler equations has the homogeneity property 



^(0= c(0e. 



(16) 
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where Cj the speed of sound in the fluid. This matrix admits a representation of the form 



C{Q)=RhR- 



(18) 



with 
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We decompose A into positive and negative parts, by defining 

v"^ = max(v, 0) , v~ = min(v, 0) , 
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and we obtain 



A^ = 






































v+ 

















v+ 




















Vr 



,A- = 



Vr 



V 










V" 












V- 





(21) 



Note that in the last expression we have assumed a low Mach number for the fluid ( |v| « Cg ), 
as expected in atmospheric dynamics. Thus, the intermediate state is defined as 
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where 
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The same argument can be appHed at the boundaries; the tests performed in section [4] considered 
soHd wall boundary conditions, which were implemented by defining ghost cells values 2o and 
23 where, 



Po 



pi, uo = ui,vo = -vi , Wo = wi , 6^0 = 



P3 = P2, U3 =U2,V3 = -V2 , W3 = W2 , 63 = 02- 

This final definition allows us to write a direct expression for the 2.5D model: 
dtQi + d,r(Qi)+ d,^(Qi) = 51(21,22), 

^.22+ ^^^(22)+ ^z^(22) = <S 2(21,22), 



(25) 
(26) 

(27) 
(28) 



<Si(2i,22) = Aj{^^(2i)+ ^"(22)- ^^(2o)- ^"(2i)}+ <S(2i), (29) 

<S2(2i,22) = A);{^^(22)+ ^"(23)- ^^(20- ^"(22)}+ <S(22). (30) 

As we previously pointed out, the resulting dimensionally-reduced model maintains the same 
structure for the physical fluxes in the x and z-directions, but has transformed the flux Q into a 
"reactive" source term that generates a coupling between the layers. The next section deals with 
the numerical approximation of this coupled system of balance laws. 



3. The numerical scheme 



In this section we present a finite volume scheme for system of balance laws derived in the 
previous section. For the sake of simplicity, we develop the scheme for a single set of equations 
of the form, 

dtQ^ d,r(Q)^ d,9{(Q)= 5(2), (31) 



being eqns. ( 27 )-( 28 ) a particular case with the augmented state Q = [Qi , 22]^- We first indicate 
that our strategy will be based in a splitting scheme, as it suggested in ||9J given the flux choice 
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that we will make. Thus, we will first establish a numerical scheme for the system of conservation 
laws, 

dtQ+ d,r(Q)^ d,9{(Q) = 0, (32) 
to be combined with a procedure for the resolution of the source term dynamics 

dtQ= S(Q). (33) 

In order to approximate eq. ( [32] ), we begin by meshing the spatial domain flx,z ii^^ uniform 
control volumes flij = [x/_i/2, Xi+1/2] x [Zj-1/2, Zj+1/2] of size AxAz. Inside every control volume 
we average with respect to x and z leading to the semi-discrete scheme 

dQi j(t) 1 1 

= -Yx^^i-mj - Fi-mj) - - H,j.,/2) ^ hj(Q), (34) 

where, 

I I rxMii rzj+i/i 
QiJ=irir\ Q(x,z,t)dzdx, (35) 

FM/2,j = ^ r ^ iQixM/2, z, 0) dz , (36) 

I rxi+i/2 

Hij^i/2 = T- *H (Q(x, Z;>i/2, 0) dx. (37) 



We approximate the expressions in eqns. d36])-(37 ) by suitable Gaussian quadrature formulas, 



Fi+i/2,j ~ ^ X ^ (Q(^i+i/2,ZGp, 0) dz , (38) 

ZGp 

Zj+i/2J))dz, (39) 



2 • 



where xcp and zg/? are prescribed Gauss points with corresponding weights Wxfj^ and w^^^ re 



spectively. The computation of eqns. (38)-(39) is performed via a high-resolution approach 



that makes use of a WENO reconstruction procedure; after this step is completed, a polyno- 
mial of prescribed order is obtained at every cell, and therefore, at every cell interface, accurate 
flux calculations can be performed by taking extrapolated boundary values. We briefly describe 
the WENO reconstruction procedure that is used in this article; we opted for the technique de- 
scribed in (TOl [m in its third order (quadratic reconstruction) version. This technique makes 
extensive use of the structure of the reconstruction procedure in one dimension, adding some 
additional mixed terms ("cross terms") that are efficiently computed by reduced stencils. It is 
an optimal and easy way to implement the algorithm for achieving high-order reconstructions 
in 2 and 3 dimensions; it also defines an unique polynomial in every cell, which is particu- 
larly useful when space dependent source terms such as viscosity are considered. At a given 
time t (the subscript indicating time is omitted throughout this derivation), given the set of av- 
eraged values {Qij} for the whole domain, at every cell, the reconstruction procedure seeks a 
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quadratic expansion upon a linear combination of Legendre polynomials rescaled in local coor- 
dinates (x,z) = [-1/2, 1/2] X [-1/2, 1/2] expressed in the form 

e(x, z) = go + QxPi(x) + QxxPlix) + Q,Pi(z) + QzzPliz) + Q,,Pi(x)Pi(zX (40) 



Except for the last term in ( [4Q| , every coefficient can be computed by performing a dimension-by- 
dimension reconstruction, which we now illustrate. We assign the subscript "0" to the cell where 
we are computing the coefficients, other values indicating location and direction with respect 
to Qo (note that the notation is coherent with the fact that the first coefficient in the expansion 
Qo, holds Qo = Qij, i.e., the centered value). Next, for this particular problem we define three 
stencils 

S ' = {e_2, Q-u Qo) , = {Q.u Go, Gi} , = {Go, Gi, G2} , (42) 
and in every stencil we compute a polynomial of the form 

Q^Hx) = Qf + QfPi(x) + Qf,P2(x) / = 1, 2, 3. (43) 

The coefficients are given by 

S' : Gi'^ = -2G-1 + G-2/2 + 3Go/2, Gi'i = (G-2-2G-i + Go)/2, (44) 



: Gf = (Gi - G-i)/2, GS = (G-i-2Go + Gi)/2, (45) 

S' : G?^ = -3Go/2 + 2Gi-G2/2, G?i = (Go - 2G-i + G2)/2. (46) 
For every polynomial we calculate a smoothness indicator defined as 

/5« = (e«)%H(2«)\ (47) 
leading to the following WENO weights: 

Y.l,a(^ (6 + /5W)'' 

where 6 is a parameter introduced in order to avoid division by zero; usually e = 10"^^. The 
scheme is rather insensitive to the parameter r, which we set r = 5. The parameter A is usually 
computed in an optimal way to increase the accuracy of the reconstruction at certain points; we 
opt for a centered approach instead, thus /l^^^ = /l^^^ = 1, while /l^^^ = 100. The ID reconstructed 
polynomial is given by 

Q(x) = J^^Q^^\x) oj^^^Q^^\x) J^^Q^^\x). (49) 

Next, a ID reconstruction in the z direction is performed in a totally analogous way. Finally, we 
address the computation of the mixed term Gxz, which is calculated in a 2D fashion. Keeping the 
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same convention regarding location subscripts as in ID, fTOl considers 4 formulas for the cross 
term upon taking all the moments around the cell. The expressions for the cross term are: 



Qi'^ = Qi,i-Qo,o-Qx-Qz-Qxx-Qzz. (50) 
eg> = -Gi,-i + 20,0 + - Gz + Qxx + Qzz, (51) 
-G-1,1 + Go,o - Gx + Gz + Qxx + Gzz, (52) 



-1(3) 



-1(4) 



G-1,-1 - Go,o + Gx + Gz - Qxx - Qzz^ (53) 

and the corresponding smoothness indicators are given by 

/5«=4(e«)%4(Q«f + (e»f. (54) 



Note that in the first part of the reconstruction, when the weights were computed, a larger subop- 
timal weight was assigned to the central stencil, which is a way to ensure stability and robustness 
of the algorithm by sacrificing additional order in the approximation (for more details, see lfT2ll ). 
However, for this term, the numerators assigned to the corresponding a's remains the same for 
every expression. The computation of this term concludes the reconstruction procedure, and 
now we have at our disposal one polynomial per cell that can be used to calculate values at the 
boundaries or inside the cell. The next step in our numerical scheme consists of the calculation of 



the numerical fluxes ([38])- ([39]), which will use extrapolated boundary values of the reconstructed 
polynomials. Rather than the use of the classical WENO scheme (as in |[T3]| or |[T4l ). which 
performs this calculation via a first order flux, we opt for the WENO-TVD approach described 
in [6]. We make use of the 2D extension of the flux-limiter-centered scheme (FLIC) approach 
presented in (TKTSl, which is a second-order, centered and non-oscillatory flux. In our case, it 
consists of a flux-limited version of a generalized Lax Wendrofl' flux, using as a low-order flux 
the GFORCE (generalized first order centered) flux 1 16], which can be interpreted as a convex 
combination of Lax-Friedrichs and Lax-Wendrofl'-type of fluxes: 

where 

= ^^f' {QUj^ ef.,,,) = o.F^ru2j - (1 - ' (56) 

F^nj = \[r {QLn,) - r (ef,,,,,) - \ ^ (ef,,,, - ef,,,,)) , (57) 

^^,%2j = r {Qlu2,j) ' (58) 

ai.u2j = \ {Qlu2j + ef.1/2,) - ^ ( ^ (ef.1/2,) - r (Gf^.^^J) ■ (59) 

The parameter oj varies between and 1, and is chosen in a compatible manner with the CFL 
number in order to ensure monotonicity; throughout the numerical experiments presented in this 



article, we will restrict ourselves to the usual FORCE flux, i.e., co = 0.5. We have omitted the 
formulas for the remaining cell boundaries, but they can be derived in a straightforward manner. 
Also note that even though the formulas are written along the boundary ' / + 1 / 2, / , the use of the 
Gaussian quadrature formula will replace the axes'/ by Gauss points and therefore this subscript 
must be understood in that sense. It is important to notice that so far we are deriving expressions 
for the semi-discrete approximation of the system of conservation laws, however, the fluxes 
include the parameter which arises from the averaging operators that originate these fluxes. 
Thus, in the spatial discretization of the system, the time stepping enters just as a parameter. At 
the end of the derivation of the scheme, when we present the time discretization of eq. ([34]), A^ 
will be considered as "marching parameter" in the sense that its inclusion in the formulas will 
generate an updated state in time. The function i//i+i/2,j = 1/2,7(^^1/2 7' ^1/2 P ^ limiter; 
a slight variation of the usual limiters has to be considered in this context since we use a centered 
flux instead of an upwind approach (the reader can refer to 1 15 , Ch 13.] for more details); in our 
case we mainly use the VanLeer limiter, which on its centered version reads: 

[0 if r < 0, 

^(^)= 1^ ifO<r<l, ^8 = \i^r (60) 

2r(l-0„) .r ^ 1 i Kl 

where c corresponds to the Courant number which depends on the problem. To implement the 
limiter in the "conservation paradigm", the flow parameter r is defined via the total energy of the 
system, 

e = c^Ott + -(u^ + + w^) + gz. (61) 
The left and right flow parameters are then defined as, 

— — p^ 

_ ^i-ll2,j ^i-l/2,j P _ ^i+3/2,j ^i+3/2,j 



1/2,7 ^i+y2J 1/2,7 1/2,7 



and finally, 



<A/+i/2,7 = min(j/^(rf^i/2,7)' ^(^+i/2,7))- (^3) 
The above described procedure starts with a set of averaged values and ends with a numerical ap- 



proximation of the space operators involved in eq. (32). The resulting scheme is still continuous 
in time, and we conclude this section by discretizing this operator in a manner that is consistent 
with the choices that we have made in the generation of the space discretization operator. At a 
given starting time we begin by considering the semi-discrete scheme 

dQijit) 



dt 



LijiQl (64) 



bringing the system to a final state f'^^ with a time stepping A^. In order to preserve high- 
order and non-oscillatory properties in time, we consider the well-known family of explicit TVD 
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Runge-Kutta schemes ifTTl , in particular its third order version 



e;;^ = oi^^ + Mum^^), (65) 

= !q« + lAfL,..(e"^^), (66) 

^i,j 4^^7 4^«j 4 ^'J^^ij ^' ^ ^ 

(68) 

We conclude this section with the inclusion of the source term. The source term appearing in 



eq. (33) will not depend on space nor space derivatives, and therefore it can be averaged in 
space and solved in the same manner as the above presented time discretization, by replacing 
Lij(Q^j) by S (Q^j)- If we denote by the £(A0 the fully discrete operator that brings the system 
of conservation laws ( [32|) units ahead in time, and by (5 (AO the fully discrete operator that 
updates the source term |33|) in A^ units, we preserve, at least, second order accuracy in time by 
implementing a Strang splitting ITSl in the form 

Q^f = e(At/2mAtmAt/2)Qlj. (69) 

In this way we have derived a fully discrete, high-order and total variation diminishing nu- 
merical scheme for treating a system of balance laws arising from dimensional reduction of the 
original set of 3D primitive equations. In the next section, numerical experiments are performed 
with the aim of understanding the consequences of the dimensional reduction and the coupling 
effects, and its ability to represent atmospheric phenomena in a plausible way. The test are mo- 
tivated by standard study cases for non-hydrostatic model development (see (H [T9j [20]|), and 
therefore a quahtative comparison with pubhshed model outputs is possible. 



4. Numerical experiments 

The purpose of the numerical experiments is to use well-known test cases that are used exten- 
sively to test dynamical cores of atmospheric models so that our results can easily be compared 
with results in the literature. There are quite a few such test cases published, and we have se- 
lected three. In addition we have constructed a new test case that is used to show the capabilities 
of the dimensionally-reduced model. 

4.1. Convective bubble in a neutral atmosphere 

This first test case, that has been previously studied in (191 ED El (among others), is per- 
formed by using the scales and settings as in [1201 . The domain is Q = [-10000, 10000] x 
[-10000, 10000] X [0, 10000] X [0, 1000], with Ax = Az= 125[m] and At chosen according to, 

= „ , ^""i , ;^ (70) 

max {\vel - c^|, \vel + Cs)\ 
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where the CFL number is set to 0.4, and vel = ^u^ + is the module of the x - z velocity field 
(this setting of the time step will be kept for the remaining tests). Both layers are initialized at 
rest, and we consider the following potential temperature profiles: 

^1 = ^1 + ^;, ^2 = ^2, (71) 

with, 

, flOcos(^) L<1, 1 r- 

^i=^2 = 30ora, e[ = \ ^ = ™V^' + fe- 2000)2. (72) 

10 i.o.c. zOOO 

The reference states 6i and 62 are in hydrostatic balance at their respective layers, and therefore 
the density p is initialized via the hydrostatic balance law, 

CpO^ = -g, (73) 
dz 

together with the equation of state for an ideal gas, 

Po ^ 

P=^^^^, (74) 

both evaluated at the corresponding reference state = 0. We implement solid wall boundary 
conditions along the domain. The underlying physics of this test dictate that, as the temperature 
perturbation is warmer than the background state, a buoyancy force will push the perturbation 
upwards. As it starts rising, it also starts experiencing horizontal expansion because of the same 
buoyancy eff'ect, eventually generating a mushroom- shaped cloud. This has been observed in 
many tests based on 2D sets of Euler equations. We observe a similar behavior in the layered 
model, as shown in figure [T] Figure [2] shows the evolution of the potential temperature residual 
between both layers, which preserves horizontal synmietry but decreases in magnitude; recall 
that the second layer is initialized without a perturbation and therefore the system tends to an 
equilibrium state between both layers. Velocity fields in figure [3] are consistent with results 
obtained in 2D tests performed in the references mentioned in the description of the test; vec- 
tor plots in figure Hlillustrate the deformation process to which the temperature perturbation is 
subjected. Finally7[5] exhibits conservation of energy of the total system, and also confirms the 
equilibrium reached by the layers. 
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Figure 1: Potential temperature for the convective bubble test case. Colormap of the first layer at t = 0, 120, 300 and 
600[s]. Ax = Az= 125[m], 160 x 80 elements. 
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Figure 2: Potential temperature for the convective bubble test case. Colormap of the residual (^i - 62) ait t = 120 and 
600[s]. Ax = Az= 125[m], 160 x 80 elements. 



90001 
8000 1 
7000 1 
6000 1 
N 5000 1 
4000 1 
3000 1 
2000 1 
1000 1 





^0.3 








9000 






0.25 


8000 






0.2 


7000 








6000 






0.15 


N 5000 






0.1 


4000 








3000 




I0.O5 


2000 






1000 




|o 






9 



-4000 -2000 



2000 4000 



-4000 -2000 2000 4000 
X 



13 



Figure 3: Velocity field colormaps at t = 120, 300 and 600 [s] for the first layer. Left: horizontal velocity. Right: vertical 
velocity. Ax = Az = 125[m], 160 x 80 elements. 
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Figure 4: Velocity field Vector plots at r = 0, 120, 300 and 600 [s] for the first layer. Ax = Az = 125[m], 160 x 80 
elements. 



9000 
8000 
7000 
6000 
N 5000 
4000 
3000 
2000 
1000 



-4000 -2000 



2000 4000 
X 



9000 
8000 
7000 
6000 
N 5000 
4000 
3000 
2000 
1000 



-4000 -2000 2000 4000 
X 



9000 [ 
8000 
7000 
6000 
N 5000 
4000 
3000 
2000 
1000 



-4000 -2000 



2000 4000 



9000 
8000 
7000 
6000 
N 5000 
4000 
3000 
2000 
1000 



-4000 -2000 



2000 4000 



Figure 5: Total energy of the system for the convective bubble test case. 
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4.2. Interaction between hot and cold bubbles 



This test is a variation of the test presented in 1 19], which studied the interaction between 
positive and negative potential temperature perturbations, modified to keep the same scale of our 



first test case; all the domain parameters are set as in test 4.1 Horizontal velocity u is initialized 
in both layers with a value of 20[ms~^] in order to make the test more stringent. A periodic 
boundary condition is set in the lateral x-direction while solid boundary walls are kept at the 
bottom and the top of the domain. Vertical velocity in both layers is set to zero. Thermodynamic 
variables are computed in exactly the same way as in the first test, but we add a cold perturbation 
to the second layer. 




As we have previously described in the first test case, in a neutral atmosphere, a warm tem- 
perature perturbation is expected to raise; by the same principles, a cold perturbation is expected 
to fall. Figure [6] shows the initial temperature for both layers. Figure |7] illustrates the experi- 
ment: placed in the same vertical axis, the warm bubble raises while the cold bubble falls, and at 
certain instant both bubbles collide starting an eddy interaction which is governed by the same 
buoyancy eff'ects. In this experiment we also included horizontal velocity and we can observe 
that the perturbations are horizontally translated with a proper speed (both layers were initial- 
ized at the same velocity) and preserving symmetry all along the experiment. Figure [8] shows a 
decrease in the magnitude of the potential temperature residual and figures 10 and [9] exhibit the 
associated velocity field, where we can observe the formation of eddies. Figure [TTJillustrates the 
evolution of potential temperature extreme values in every layer, and it can be seen that for both 
maximum and minimum, the layers have a tendency to reach equilibrium values. Total energy is 
preserved as shown in figure [12] in a similar manner as in the first test. 
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Figure 7: Potential temperature for the hot and cold bubbles test case. Colormap of the first layer at t = 120, 300, 600 
and 1000[s]. ^x = ^z= 125[m], 160 x 80 elements. 
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Figure 9: Velocity field vector plots at t = 120, 300, 600 and 1000[s] for the second layer. Ax = Az = 125[m], 160 x 80 
elements. 
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Figure 10: Velocity field colormaps at f = 300[s] for the second layer. Left: horizontal velocity. Right: vertical velocity. 
Ax = Az = 125[m], 160 x 80 elements. 



9000 
8000 
7000 
6000 
N 5000 
4000 
3000 
2000 
1000 



I 



I 



i 



9000 
8000 
7000 
6000 
N 5000 
4000 
3000 
2000 
1000 



tit 
f 



18 



Figure 11: Potential temperature extrema evolution for the hot and cold bubbles test case. 




4.3. Adjustment of shear between two layers 

This third test case aims to study the level of interaction and adjustment between the two 
layers in the model. We again consider the same domain parameters, except for a final simulation 
time of 300[5']. However, variables are initialized in a diff'erent way. First of all, horizontal 
velocity v is set to a value of 10[m5'"^] in the first layer, while in the second layer it has a value 
of -\Q[ms~^]. The horizontal velocity u is initialized in the first layer with a logarithmic profile 



«i=50(log(^ + l)) 
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(76) 



and it is set to zero in the second layer. Vertical velocity in both layers is set to zero. Thermody- 
namic reference states are different: no potential temperature perturbation is considered in any 
case, the first layer having the same neutral atmosphere as in the previous test cases, while the 
second is considered to be a stable atmosphere with 



= ^oexpl — -d, 



(77) 



where 60 = 61= 300[^] and N is the Brunt- Vaisala frequency with a value of TV = 0.01 [^"^]. 
In both cases p is computed as in the previous tests. Boundary conditions are set again to be 
periodic in the lateral x-direction and solid wall boundary conditions in z- 



Figure 13: Temporal evolution of the horizontal velocity in the x-direction u, vertical section at x = 0, in layer 1 (left) 
and layer 2(right). 



N 20 
O 




o 
o 

o 25 
c 

o 20 
o 




6000 8000 10000 



2000 4000 6000 8000 



In figure [13] we present the velocity profile for the velocity component in the x-direction. We 
see that the profiles are adjusted to some logarithmic profile from above in layer 1, and from 
below in layer 2. This is as expected. The same adjustment is seen for the potential temperature 



in figure 14 the adjusted atmosphere becomes stable with a stability which is a "mean value" 
between the two layers. 
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Figure 14: Temporal evolution of the potential temperature 6, vertical section at x = 0, in layer 1 (left) and layer 2(right). 
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In figure 15 we present the time evolution of the velocity component in the y-direction, in a 
cross section at x = 0. So here we expect to see a decay of an initial shear, produced by the initial 
shear in the x-direction. And this is what figure p?] shows, a decay from a positive value in layer 
1, and from a negative value in layer 2. 



Figure 15: Temporal evolution of the horizontal velocity in the };-direction v, vertical section at x = 0, in layer 1 (left) 
and layer 2(right). 
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Note the interesting "undershoot" in layer 1, and the "overshoot" in layer 2: The profiles 
maintains their shape, before the final decay towards zero. Recall that this test is without explicit 
friction or parameterizations of the shear layer, so it is basically thermodynamics (6 ^ P ^ u) 



that controls the process. That the adjustment works well is also illustrated in figure 16 which 
shows the conservation of total energy and the adjustment of the energy in the two layers. 



Figure 16: Total energy of the system for the wind shear test case. 
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This test case can be extended to have a rising bubble in one of the layers as well as other per- 
turbations to simulate front behavior; one can also introduce parameterizations of the shear layer, 
for example by considering a simple turbulence model (a Smagorinsky model is one choice). 

4.4. Inertia- gravity waves through a periodic channel 

The last case that we study is a variation of the inertia-gravity waves test case proposed by (H. 
The purpose of this test is twofold: first, we qualitatively compare the behavior of our model for 
the generation of wave motion in a nonhydrostatic scale, but we also study the nonlinear eff'ects 
of the coupling. We set Q = [0, 300000] x [-10000, 10000] x [0, 10000] x [0, 3000], with 
Ax = Az = 500[m]. In a first model run, the horizontal velocity u is set to a value of 20[m^"^] 
for the first layer, while any other velocity variable is set to zero. Both layers are initialized with 
the same stable atmosphere profile as in the previous example, and a temperature perturbation is 
added to the first layer, in the form. 




with Op = 10[^], Xp = 100, 000[m], and ap = 5, 000[m]. In a second model run, we leave the first 
layer at rest and add a velocity u = -20[ms~^] to the second layer. We remove the temperature 
perturbation from the first layer and assign it to the second, only changing the parameter Xp = 
200, 000, thus expecting a similar behavior as in the first model run but in the opposite direction. 
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A third and final run considers a combination of the previous model runs, i.e., one perturbation 
in every layer, the first layer with a positive horizontal displacement and the second layer with 



a negative horizontal displacement. Figure 18 shows the results of the first and second model 



runs, initial condition and result at ^ = 2500, first run in the two upper panels and second run in 
the two lower panels. These results are in very good agreement with the results in |4 | and other 



test results in the literature. In figure 19 we present the results of the third model run. We see 



immediately that these results are not superpositions of the results in figure 18 so we see the 
eff'ects of a nonlinear interaction between the two inertia- gravity waves. In figure 20 we present 



the velocity fields for the nonlinear interaction case. We see time-dependent fan-like structures 



typical of nonlinear wave interaction and almost perfect symmetry of the fields. Figure 17 shows 
the evolution of extremal values for the horizontal velocity in every layer. It can be seen that both 
layers tend to an equilibrium, while an interaction between maximal and minimal values occurs. 
It is not easy to explain all the details in the results, for example the development from t = 2000 
to t = 2500, since the interaction is purely nonlinear. However, it is possible to see eff'ects of 
focussing, which is a purely nonlinear wave phenomenon. In the solutions the symmetry is well 
preserved indicating that the is performing in a good manner; conservation of energy is shown in 
figure|2T]. 



Figure 17: Horizontal velocity u extrema evolution. Inertia-gravity waves case; combined test run. 
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Figure 18: Inertia- gravity waves. Potential temperature colormaps at t 
500[m], 600 X 20 elements. 
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Figure 19: Inertia-gravity waves. Potential temperature colormaps for the first layer at f = 0, 1000, 2000 and 2500[5] in 
the combined test. Ax = Az = 500[m], 600 x 20 elements. 




25 



Figure 20: Inertia- gravity waves. X - Z vector field plot for the first layer at f = 0, 1000, 2000 and 2500 [5] in 
combined test. Ax = Az = 500[m], 600 x 20 elements. 
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Figure 21: Total energy of the system for the inertia-gravity waves case; combined test run. 




5. Summary and concluding remarks 

We have derived a set of dimensionally-reduced fully nonlinear primitive equations to model 
atmospheric dynamics. Our goal was to show that the 2.5D model is able to reproduce well- 
known atmospheric phenomena properly. A DG approach was used for the dimensional reduc- 
tion, whereas a WENO-TVD method was used for discretization. We think that finite volume 
methods are the most promising discretization method for atmospheric models since it ensures 
discrete conservation. If one uses, as we have done, equations in conservative form with con- 
served variables, the discretized equations should model the conserved quantities in the atmo- 
sphere with high accuracy. The numerical experiments with extensively used test cases for at- 
mospheric model cores have shown that the dicretization is accurate, stable, energy-conserving 
and do not produce spurious oscillations. No artificial diff'usion or any other form of smoothing 
is used. We think that this is a strong indication that our set of equations are well-posed, and 
that a higher order discretization is indeed valuable. Our approach can be extended in several 
directions: we could use higher order elements in the dimensional reduction and/or even higher 
order finite volume methods. The purpose of this would be to investigate the accuracy versus the 
computational costs. This can easily be done for well-known test cases where we have results 
to compare with, but we would also like to use more complicated tests cases from real weather 
situations. The latter test cases could be used to investigate how one can model certain weather 
phenomena where one, from a physical viewpoint, have a certain insight in the behavior. There 
is also an option of using other methods than the DG method for the dimensional reduction, but 
our opinion is that this method is good natural choice, because it poses no difficulty for handling 
nonlinear equations. Also on the computational side, the actual computational costs has not been 
explicitly addressed in this paper. However, it is worth to mention that every block of the code 
can be completely vectorized, while the dimensionally-reduced equations have a natural form of 
parallelism by parallelizing over the layers. This is an important topic since the existing dynami- 
cal cores of atmospheric models which we would like to compare with in term of computational 
cost, are parallelized. 
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